Nonequilibrium patterns and shape fluctuations in reactive membranes 
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A simple kinetic model of a two-component deformable and reactive bilayer is presented. The two 
diflterently shaped components are interconverted by a nonequilibrium reaction, and a phenomeno- 
^ \ logical coupling between local composition and curvature is proposed. When the two components 

O ■ are not miscible, linear stability analysis predicts, and numerical simulations show, the formation of 

' stationary nonequilibrium composition/curvature patterns whose typical size is determined by the 

reactive process. For miscible components, a linearization of the dynamic equations is performed in 
order to evaluate the correlation function for shape fluctuations from which the behavior of these 
systems in micropipet aspiration experiments can be predicted. 



PACS numbers: 87.16.Dg, 07.05. Tp, 82.45.Mp 



C3 . I. INTRODUCTION 



Polar lipids self-assemble and orient, with the hydrophilic portions facing water. The water may be sandwiched by 
two lipid layers, as in the black spots in soap bubbles that develop at locations where the two layers are closer than 
the shortest wavelength of the visible spectrum and that occur just before the bubble bursts, or the water may be 
^ outside of the two lipid layers, as in biomembranes. While it is known that the structure of a biological membrane 
is far more complex than a simple lipid bilayer because of embedded proteins, cholesterol molecules, and ions, 

to name just a few of the components that provide the full functionality to this structure that is crucial to the life 
J> ! of t;he cell, it is nevertheless acknowledged that the lipid bilayer is the basic structural unit of all cell membranes. 
' An understanding of this simpler system is therefore extremely important in an effort to shed light on the properties 
' and functioning of active transport, signaling, and adhesion in cells. Furthermore, artificial lipid bilayers are widely 
. used in a number of nanotechnological applications ranging from solar energy transduction and biosensors to drug 
development. 

A lipid bilayer is highly flexible and liquid-like (as is a real membrane). It can therefore not be viewed as a static 
inert boundary but must be recognized as a dynamical structure In the case of a biological membrane, lipid 
bilayers serve as quasi-two-dimensional solvents for proteins and all other components, and are intimately involved in 
■ many biochemical processes. Moreover, its ability to change its shape (curvature) is an essential property of the lipid 
bilayer since the formation of vesicles and the permeability properties of the cell depend on it. 

From the modeling viewpoint, it became recognized during the 90's that some internal degrees of freedom are 
' necessary to understand the large variety of conformational changes found in cell and synthetic membranes. In 
Ch particular, the local composition of the bilayer can crucially affect its local curvature, and some models were developed 
O on the basis of this idea ^^^^J^^^- However, these approaches considered membranes as equilibrium systems. 
In a biological context, this hypothesis is at best hopeful. 

Nonequilibrium conditions are ubiquitous in nature, and for this reason the out-of-equilibrium behavior of mem- 
branes, both in cellular systems and in laboratory-prepared vesicles, has attracted much attention over the past few 
years. The first successful approach to the study of nonequilibrium membranes was introduced by J. Prost et al. [l0| . 
cd ' Roughly summarizing their modeling approach, they suggested that some externally activated components (intramem- 
brane proteins) act as "pumps," generating forces on the membrane that locally change its curvature. Variations, 
improvements and sequels of the initial model [llL ITsL as well as related experimental studies |0, 0| to a 
large extent complete the understanding of the nonequilibrium behavior of bilayers with inserted active components. 

Our interest in this work lies in a different nonequilibrium that may arise from externally activated chemical pro- 
cesses involving the transformation of the elementary membrane lipid components. As an example, nice experiments 
by P. G. Petrov et al. [3 

show the effect of an ongoing photocontroled chemical reaction on the curvature of synthetic 
giant vesicles. Other evidence of chemically-induced shape transformations is found in the nervous synaptic process: 
a reaction that interconverts two differently shaped constituent phospholipids of the membrane plays a crucial role in 
the fission of vesicles in nervous cells 17, 18] . Recent experiments also point to the importance of high-curvature lipids 
in many cell processes such as membrane fusion [Tgj . So far, no models for this kind of nonequilibrium situation can 
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be found in the literature. Our aim is to gain physical insight into the role of a generic nonequilibrium reaction acting 
on a membrane that is, in turn, described through a curvature/composition coupling. We show how the reactive 
process endows the membrane with a rich pattern formation phenomenology and specific characteristics of its shape 
fluctuations. 

We present a simple kinetic model of a two-component deformable bilayer with a chemical reaction. We model 
nonactive membranes where the nonequilibrium nature of the system is caused by a chemical reaction that intcr- 
converts two differently shaped components of the bilayer rather than by inserted active components. In the model, 
the membrane is composed of two species: lipid A, which is assumed to be coned-shaped, and lipid i3, which has an 
inverted cone shape (see Fig.^. We consider the simplest scenario where the outer layer of the membrane is composed 
of A and B lipids, whereas the inner layer is composed of a single component without any curvature effect. We also 
prescribe that lipids do not move between the inner and outer layers. Both of these simplifications are made in most 
membrane models. According to these assumptions, the membrane is simply modeled as a laterally heterogeneous 
elastic surface with an internal composition order parameter locally coupled to the curvature. 




FIG. 1: The two membrane components A (left) and B (right) have opposite cone shapes. We arbitrarily choose a positive 
curvature for A and a negative one for B. A reaction interconverts both species with k+ and k- being the forward and backward 
reaction rates, respectively. 



Within this framework, we are basically interested in two situations. On the one hand, for immiscible components, 
we show that phase separation in the membrane leads to the spontaneous development of structures involving hetero- 
geneous distributions of both composition and curvature that finally result in stationary finite-sized nonequilibrium 
domains. This instability is studied by performing a linear stability analysis of the kinetic equations, and the pattern- 
selection role of the reactive process is established. Corresponding numerical simulations will be shown to support 
these predictions. On the other hand, for miscible components we derive the correlation function for the membrane 
shape fiuctuations. The fluctuation spectrum is mainly determined by the bending energy at small surface tension, 
and here we find the same expression for nonreactive and reactive cases, but with different bending rigidities. This 
change in the behavior of the fiuctuation spectrum between equilibrium and nonequilibrium states might be easily 
controlled and observed in micropipet aspiration experiments. 

This paper is organized as follows. In Sec. ^1 the study of membranes of immiscible components is addressed. 
We propose a free energy functional, derive the kinetic equations, and perform the linear stability analysis of these 
equations. Numerical simulations are carried out and some representative results are shown. In Sec. IIIII bilayers 
of miscible components are analyzed and the change of the height fluctuation spectrum between the reactive and 
nonreactive situations is established. We conclude with a brief summary in Sec. IIVI 

II. MEMBRANES OF IMMISCIBLE COMPONENTS 
A. Model and analytical results 

The membrane is defined as a two dimensional surface with a concentration difference order parameter and a local 
extrinsic curvature H. The rigidity of the membrane leads to an elastic energy contribution ^ J [H — Hsp{(f>))^dxdy 
to the total energy, where k is the bending rigidity modulus, and the spontaneous (equilibrium) curvature Hsp((j)) 
refiects the shape asymmetry between the two lipid components. For simplicity, we adopt a linear dependence on 
(j), Hsp = ipHf), with Hq > according to the schematic in Fig.^ In the Monge parametrization 0| a deformable 
surface is described by [x, y, h{x, y)), where h{x, y) is the displacement (height) field for the local separation from the 
fiat conformation. This representation is valid for surfaces that are nearly fiat with only gradual variations of h, and 
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allows the approximation H h- As a function of these variables, the proposed energy functional reads 



dxdy, (1) 



2^ V V ^' 2 

where the first three terms correspond to the typical Ginzburg-Landau expansion responsible for phase separation (a, 
/?, 7 > 0), with an equilibrium concentration difference = ±-\/q;//3, and a typical interface length C, = \p~ffoL. For 
self-assembled free membranes, the surface tension contribution (||V/ip) in the free energy can be neglected, and we 
have not included it in Eq. 

The kinetics of (/) follows a conserved scheme |23| plus the reaction contributions, 



r(</'-0o), (2) 



where F = fc-)_ + fc_ and 0o (fc- ^ + fc-)- fc+ ^.nd fc_ are the forward and backward reaction rate constants, 

respectively. Considering a permeable membrane (i.e., ignoring hydrodynamic interactions) we adopt the following 
relaxational equation for the evolution of the height field, 

where A is a mobility parameter proportional to the inverse of the typical relaxation time r^. 

The kinetic equations can readily be adimensionalized: energy is measured in units of ksT, time in units of r/j, and 
length in units of \/Dt\. In terms of the new dimcnsionless parameters, the kinetic equations become 

^ = -a)v20 + 3/3(/)2v2</. + 6/3(?;.|V0|2 

- 7VV - kHq^^K - F - (/-o) , (4) 

dh . o 

— = -kV*/i + KiJoV^^. 
ot 

At thermal equilibrium Eqs. Q describe the spinodal decomposition of two immiscible components. Due to the 
composition/curvature coupling Hq, both fields will form complementary patterns. As phase segregation progresses, 
membrane regions with positive (negative) (j) deform in such a way that the curvature become positive (negative), as 
shown in Fig. |21 In the absence of reaction, this coarsening process does not end until there is complete segregation 
into two large domains. A nonequilibrium reaction such as the one we propose converts one species into the other. 
This amounts to a large-scale mixing mechanism that counteracts the short-scale ordering effect of phase separation. 
Therefore, the segregated structures grow only until mixing and ordering effects compensate, resulting in a stationary 
pattern. These types of nonequilibrium patterns are also found in other systems such as polymer blends as 
well as in monomolecular adsorbtion on metal surfaces |3i US j and have to be distinguished from typical Turing 
patterns |2^. Even though they emerge from the same kind of instablity (see below), Turing patterns are expected in 
systems with species with different diffusivities. In the model presented here, the domains result from the competition 
between a local thermodynamic affinity of equal species and a nonequilibrium reaction mixing effect. As we will see, 
linear stability analysis and numerical simulations support this idea. 

The stationary uniform state corresponds to 4> — (t>o and arbitrary h. The linear stability of these uniform solutions is 
tested by adding small plane- wave perturbations of wave numer q and linearizing Eqs. This procedure determines 
the 2x2 linearization matrix C, with the following coefficients, 



£11 = -g^ [(^^2 _ ^ ^ 3^^2^ ^ 
£12 = -nHoq^ 
C21 = -KHoq^ 
C22 = -Kf?**. 



21 



(5) 



The eigenvalues ujq of the Jacobian associated with the matrix C correspond to the linear growth rates of the 

perturbations. Solving the eigenvalue problem we obtain ojq = ^ (Tr[C] ± ^/A\£]^ , where A[C] = Tr[C]^ - 4:Det[C]. 

At the instability boundary, Re{ijjq) vanishes for one finite wave number that is defined as the first unstable mode. 
If the imaginary part of LDq is not zero at this wave number, we have a wave bifurcation. The condition for this 
bifurcation is obtained by requiring T'r[£] = and A[£] < 0. For positive k, F, a, /3 and 7, these conditions do not 
apply at any real wave number, and consequently, wave instability is not found in this model. 
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FIG. 2; Schematic representation of the composition/curvature couphng effect in an unstable membrane. 



On the other hand, if the imaginary part of the growth rate is zero at the bifurcation point, we have a Turing-hke 
bifurcation. The condition for this bifurcation is Det[C] = 0, whose analytical expression is easily obtained, yielding 
the following condition for the model parameters. 



A phase diagram is shown in Fig. |21 In the unstable region, Turing-like patterns emerging from the competition 
between phase separation and reaction are predicted. In the stable phase, although the two components are not 
inherently miscible, the reaction completely mixes the system, and the bilayer becomes stable and essentially flat. 

In the unstable phase, F < Fc, there is a range of unstable modes. Comparing with the equilibrium case (F = 0) 
where a continuous range of modes starting at g = is unstable, reaction stabilizes the long wavelength modes so 
that, as we anticipate, the long time distribution of the system results in a stationary nonequilibrium pattern with a 
finite size (see further discussion in Sec. IIIBII . Notice in Fig.0]that increasing F reduces the range of unstable modes 
progressively until F reaches the marginal value in Eq. (jSJ , above which there is no longer any instability. The limits 
of the range of unstable modes are given by 



which are independent of curvature parameters. Curvature reduces the unstable mode growth rates (see line with 
circles in Fig.Q but without changing either q± or the marginal condition ©. The effect of curvature is exclusively 
related to the kinetics of the phase separation process (see below) . 

Before numerically solving the kinetic equation, we can seed light on the expected pattern formation process by 
performing a weakly nonlinear analysis using the amplitude equations technique. This analysis allows us to compute a 
solution for our problem near the bifurcation threshold, find the universality class of the pattern formation mechanism, 
and explain some properties such as the spatial arrangement of the patterns found in the numerical simulations. 

For simplicity, we restrict the calculation of the amplitude equations to one spatial dimension. In spite of this 
restriction, we will be able to predict the kind of spatial arrangement found in two dimensions by using the universality 




(6) 



and for the wave vector of the first unstable mode. 




(7) 




(8) 
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FIG. 4: Dispersion relation functions ujq at different values of V. The other parameters are held fixed at q = 7 = 1, (/jq = and 
K = Ho — 0, except for the curve with circle symbols (labeled with '& curv' in the legends) that corresponds to k = 10 and 
Ho = 0.2. 



properties of the amplitude equations. Details of the derivation of the amplitude equations are presented in the 
appendix. The analysis reveals that a solution for Eqs. Q near the bifurcation to pattern formation reads 

(j) {x, i) = 00 + [^exp (iqcx) + c.c] , 
h{x,t) = h + [B exp {iqcx) + c.c] , 

where c.c. stands for the complex conjugate, and the amplitudes A and B satisfy the equations 



dtA = (r, ~r)A 



(« - ml 



2^2 



3..A^'J^^^^±^^9..B, 



dtB = Hand.^A + 



(9) 
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Moreover, since B is purely relaxational, we can adiabatically eliminate it and reduce ® to a single evolution equation 
for the amplitude A, 



Therefore, the amplitudes satisfy the real Ginzburg-Landau equation. If F > Fc then the amplitudes relax to zero 
and a homogeneous state (0 = 0o and arbitrary h = h) is obtained. However, if F < Fc then the amplitudes 
reach a stationary value and a pattern develops if qc G M. Notice that the amplitudes present a negative-positive 
aspect, that is, at sites where A reaches a maximum (minimum) B reaches a minimum (maximum). This is simply 
a consequence of the particular selection of sign for the spontaneous curvature Hq. With regard to the nonlinear 
term, ^ note first that its coefhcient is always negative. As a result, the nonlinearity always plays a stabilizing 
role, the bifurcation to pattern formation is supercritical for all values of the parameters, and no hysteresis may 
occur. Secondly, the nonlinear term provides information about the relevant modal interactions. At this point, we 
can take advantage of the universality properties of the amplitude equations to infer the spatial arrangement of the 
pattern in two dimensions. Since the Swift-Hohenberg model also shares the same universality class when performing 
an amplitude equation analysis, namely, it also reduces to the real Ginzburg-Landau equation psL l2^. the pattern 
formation mechanism is the same in both models. It is well-known in that case that the modal interactions in two 
dimensions are such that if the inversion symmetry is preserved then roll-like patterns develop. On the other hand, 
if that symmetry is not fulfilled, a hexagonal structure appears. Note that in Eqs. the inversion symmetry, 
{(f) — > —(/), h — > —h} is only satisfied if 0o = 0. Thus we expect a roll-like pattern in that case, while hexagons will 
develop otherwise. 



Numerical integration of Eqs. has been performed in two dimensions using an explicit Euler scheme in a square 
lattice with periodic boundary condition. Small random perturbations around (j) = and /i = are implemented 
as initial conditions. The coordinate step Aa; was chosen equal to 1, and the time step was usually Ai — 10^"* to 
assure good numerical accuracy (length and time in dimensionless simulation units). The numerical results presented 
in this section correspond to highly immiscible components (deep quench, a = \ and 13=1, leading to an equilibrium 
value of (f)^q = and an interface thickness of the order of the space discretization {C, = 1, which leads to 7 = 1). 
The bending rigidity modulus is taken equal to 10 (in units of ksT) p^ . and we consider two constituent lipids of 
very different shapes by setting Hq = 0.2. All our numerical results are consistent with the predictions of the linear 
stability analysis and the amplitude equations. Representative numerical results are presented below. 

In Fig. 13 we show the simulation results for three different situations. The first row corresponds to small F = 0.05 
and a critical quench {(po = 0), showing the development of a laberynthinc pattern that is still evolving at t = 2000 
(the longest time for the snapshots shown in the figure). The coarsening process, however, is progressively slowed 
down later on. The ^-field profiles of these domains (horizontal cross-sectional cuts through the pattern) reveal regions 
where 4>eq is -1-1 or —1 connected by abrupt boundaries that indicate the short spatial range over which the equilibrium 
order parameter changes from one value to the other. This is a signature of the fact that the phase separation process 
is dominant for this value of F. The second and third sets of panels correspond to large F = 0.2 (the marginal condition 
for the selected parameters is Fc = 0.25) showing laberynthinc (critical quench in the second row) and quasi-hexatic 
droplet-like (off critical quench, c/iq = 0.1, in the third set) patterns. A spatial Fourier transform of this stationary 
pattern shows a clear hexagonal structure. These domains are already stationary at times longer than t = 2000, and 
can indeed be considered as nonequilibrium stable phases of the system, involving both composition and curvature 
modulations. Their 0- and /i-field profiles have a smooth harmonic shape due to the strength of the reactive process. 
In both cases, since we are close to the bifurcation boundary, small local deviations from the stationary values are 
obtained. We especially monitor the variations of the height field and find that (|V/i|) < 0.05 is satisfied, so that the 
model has a real physical correspondence. 

In order to assess the kinetic ordering process more quantitatively, we monitor the domain size L(t) computed from 
the composition correlation functions (0(r', t)(/)(r' -I- r,t)) (the same results are obtained using the height correlation 
functions). Here, the brackets (• • • ) indicate not only an average over orientations of r and over surface positions 
r', but also over different realizations of random initial conditions. As has been reported in other studies, curvature 
considerably slows down the coarsening segregation process. This is specially evident when reaction is absent, and for 
this situation different kinetic approaches have been invoked. A mean-field kinetic scheme by Taniguchi leads to 




(10) 



B. Numerical results 



FIG. 5: Three sets of selected patterns resulting from the numerical simulation based on Eqs. First row: F = 0.05, (j>o = 0. 
Second row: F = 0.2, (f)o = 0. Third row: F = 0.2, <j)o =0.1 (off critical quench), showing an array of droplet-like domains 
rich in the minority species. The first three panels in each row correspond to the ^-field distributions of a 128 x 128 system 
a.t t — 100, 400 and 2000 (from left to right). Darker (lighter) regions are richer in the A (B) lipid. In the last panel of each 
row the height field has been plotted in three dimensions for a 64 x 64 portion of the corresponding system at i = 2000. Some 
exaggeration along the vertical direction has been applied. 

extremely slow laberynthine stripe growth obeying L(t) ~ with r — 0.1, instead of the usual spinodal decomposition 
growth exponent of 1/3 [s^]. Later, Monte Carlo simulations [11,0] showed how, in the long time evolution, those 
stripes break up into disconnected buds that subsequently diffuse and coalesce. In our model, when the reactive term 
is removed, similar results as those presented in Ref. are obtained, and no stripe breakup is observed. The reason 
for the disagreement between the analytic and Monte Carlo schemes is probably due to the fact that continuum 
models subjected to a specific parametrization of the surface such as given in Eqs. Q do not allow for overhangs that 
are surely crucial in the late stages of membrane phase segregation. 

However, the cases in which there is no reaction or in which the reaction is weak are not of interest for us, since in 
these cases the system evolves in such away that large gradients of the height displacement field occur. As noted above, 
when this happens the Monge surface parametrization is not valid and the model, although mathematically robust, 
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does not describe the physical behavior of any real system. In the reactive cases, however, the slow kinetics is still 
observed for the times prior to stationary pattern formation. This is observed in Fig.|^ where a set of curves showing 
the domain size L{t) is presented for several values of T. The situations with and without curvature are compared for 
each case. Notice how the deformable systems evolve more slowly than the nondeformable ones, although the same 
final stationary size is achieved. There is no theory to explain such a slowing down effect, but some hand-waving 
arguments may explain the physical reasons for this behavior. Monitoring of the different energy contributions in 
Eq. JQ) indicates that at early times (when the kinetics with and without curvature are still quite similar) interfacial 
energy due to the rapid formation of composition domains is much larger than any other energy contribution. In 
other words, phase separation proceeds and membrane curvature follows the composition change. At this stage the 
reduction of the interfacial energy governs the coarsening process, leading to the 1/3 Lifshitz-Slyozov exponent. Later, 
at intermediate times, when composition domains become larger, the height order parameter is no longer able to keep 
up with the phase separation process and some curvature energy is stored, causing the whole coarsening process to 
slow down. 

1.5 
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FIG. 6: Log-log plots of L(t) vs t for different values of F. The values of L{t) have been computed as an average over 10 runs 
of the correlation functions in 128 x 128 systems. In all cases (j>o = 0. The curvature parameters are set to k = 10, -ffo = 0.2 
for the thin lines, and k = Ho = for the thicker lines. 



The final size of the nonequilibrium stationary domains, Lf — L{t oo), is determined by the reaction parameter. 
The dependence of the final pattern size on the reaction parameter, Lf ^ L^*, has been largely discussed in the 
literature |2^|^|33|. The results of our model also reproduce the two limiting behaviors: s = 1/4 for large F (close 
to its marginal value), and s = 1/3 for small F. In Fig. |7| i/ is plotted for different values of the reaction parameter. 
Linear fits for the four first and the last four data points give the slopes 0.323 ± 0.006 and 0.26 ± 0.006, respectively. 
The derivation of the exponent values for the rigid situation (nondeformable surfaces) is performed by minimizing an 
effective free energy expression (where the reactive term is included via Green's functions) in a square (small F) and a 
harmonic (large F) approximation for the stationary concentration field |3lLl32j|. In the present model, the curvature 
kinetics relaxationally follows the concentration dynamics. At longer times (when the stationary patterns are already 
achieved) no significant curvature energy is present in the system, so that we can neglect the curvature contributions 
in Eq. and the results in Refs. ^3lll3^ are recovered. 

In Eq. (PI we ignored the hydrodynamic effects due to the background fluid velocity by considering a permeable 
membrane through which the fluid drains freely as the system evolves. Solvent hydrodynamic effects [33 can be 
approached through a renormalized height mobility A = (4/i(7)~^ in the second of Eqs. ^ in Fourier space, where /i 
is the solvent kinematic viscosity. The inclusion of this g-dependent mobility does not change the main results of the 
stability analysis and numerical simulations shown so far. 

Before ending this section, we check the applicability of our results to real membrane systems. We adimensionalized 
the kinetic equations to units in which D = K = ksT = 1. Redimensionalizing them, taking the typical values 
D = 10^^ — 10^*cm^/s (for lipids in a liquid-phase membrane), /i = Icp {H2O at 20°C) and the typical unstable 
modes q ~ Qc — 0.5, we find that the size of the patterns obtained above lies within the range 2 — 20//m, which is 
accessible in giant vesicle (10 — lOO^m in diameter) experiments. 



r=0.02 

r=o.05 
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FIG. 7: Log-log plot of L/ vs. F. The slopes 1/4 and 1/3 are plotted for comparison with the large and small F regimes, 
respectively. All the points have been calculated for (po ~ 0, k = 10 and Hq = 0.2, at sufficiently long times to consider the 
growth process practically halted. We have performed an average over 10 realizations. The system sizes are 128 x 128 for large 
F, and 256 x 256 for small F. 



III. MEMBRANES OF MISCIBLE COMPONENTS: LINEARIZATION AND SHAPE FLUCTUATIONS 



In this section we consider a parameter region where the membrane is thermodinamically stable. The instabihty in 
the previous section was caused by the immiscibility of the two lipid components, while we now consider the case of 
miscible membrane constituents. Notice that the two situations could correspond to the same physical system but at 
different temperatures, below (unstable) or above (stable) the critical temperature. For the stable case we adopt the 
following free energy functional. 



dxdy, 



(11) 



where a is now negative, and the nonlinear (/3) and line tension (7) terms have been removed since they are irrelevant 
in the absence of phase segregation. An additional surface tension term has been included in order to study membranes 
under tension. 

In the equilibrium situation (no reaction), the concentration field can be integrated out (^ — 0) from Eq. Ullfl. 



leading to 



-'V^h. Inserting (f>eq in Eq. Hll|l we obtain an effective free energy for h alone, 



^eff = 



-eff 



with an effective rigidity modulus [s^ ] 



dxdy, 



(12) 



(13) 



Notice that, as a consequence of the composition/curvature coupling, the different preferred curvature of the lipid 
components acts to reduce the rigidity modulus from k to Ke//- Note also that n^ff > for negative values of a, so 
that the stability of the membrane is assured above the critical temperature. 

The dynamics for the immiscibility situation has so far been considered to be strictly deterministic. Here, however, 
in order to obtain dynamical steady state functions, we add stochastic forces to the kinetic equations and perform an 
average over an appropriate ensemble. The dimensionless kinetic equation for (j), once the reactive terms are included, 
reads 



dt 



(14) 
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where the last term corresponds to a conserving Gaussian noise with correlations (/^(r, t)/^(r', i')) = 2kBTS{r 
r')(5(t — t'). For the curvature order parameter we again adopt the relaxational equation for a permeable surface, 



dh 
'dt 



(15) 



with A = 1 in dimensionless units, and fh is a thermal equilibrium noise with correlations {fh{'r,t)fh{'r'^t')) 
2kBTS{r~r')S{t~t'). 

The kinetic equations become 



kHqW^h - r I 



dt 



(16) 



The important quantity to characterize membrane shape fluctuations is the height variance at wave number 

q, which is calculated by linearizing the kinetic equations H16() and solving them in Fourier space. If (/)(q, u) and 
h{q,uj) are the Fourier transforms of (j){r,t) and h{r,t), respectively, Eqs. H16|l can be written as 



ail ai2 

0-21 0,22 



0(q,w) 



fh 



where 



ail -- 
ai2 = -nHoq"^ 
a2i = -K.Hoq'^ 
0-22 = —nq'^ + crq^ . 

Solving these coupled equations, and using the statistical properties of the thermal noises, yields 



{h{q,uj)h*{q,uj))^ 



2kBT{aliq^ 



(flu + 022)^ + (012021 - 011022 + tjJ^ f 



On integrating over w, the height variance follows the general expression. 



(1^ 



1 



V — r\ 



A 



A' 



where v and are negative variables that read 



and 



A- 
A' 
B 
C 



2kBT {aliq^ + all) 
2kBT 

On + O22 + 2a2iai2 
(011022 - 021012)^ . 



(17) 



(18) 



(19) 



(20) 



(21) 



(22) 



In the absence of the reaction, the evaluation of Eq. (|20|l in the long-wavelength limit for a tensionless membrane 
leads to — kBT / Kef fq^ ■ This behavior is truncated and replaced by kBT jaq^ for a membrane under tension if 

only the dominant terms at small q are retained. Keeping the dominant and the first subdominant terms, one recovers 
the well known expression for the height variance in an equilibrium membrane. 
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Notice that this resuh is obtained in a much simpler way from an equilibrium average using the effective free energy 
in Eq. (P|l . 

One of the usual experiments to study equilibrium membrane shape fluctuations is based on the micropipet as- 
piration technique j33 |. In these experiments, a pressure difference is applied inside a micropipet in contact with a 
vesicle membrane. This creates a tension in the membrane that pulls the excess area AS* due to the thermal shape 
fluctuations inside the micropipet. By means of this technique, the areal strain a = AS/S is obtained experimentally 
for different values of the applied tension a. According to Eq. 123|) . the relative areal strain Aa = a — ao [ao being 
the areal strain for a reference value ctq) can be calculated analytically 36, 37]: 

A^=^ln(^). (24) 

Therefore, the slope of the logarithm of a versus Aa obtained by the micropipet technique yields a measure of the 
effective bending modulus. Note that in these experiments a certain tension cr is needed, but it has to be small since 
otherwise the Keffq^ term in Eq. H23|l would be insignificant compared to aq^. 

Now the question is, how is the shape fluctuation spectrum affected by the presence of the reaction? In other 
words, how would the nonequilibrium membranes described here behave in micropipet experiments? To answer this, 
we evaluate Eq. H20() for T ^ and keep only the dominant terms in the limit q —> 0. For tensionless membranes we 
get = ksT/ nq^ and for membranes under tension we recover = jaq^ . Keeping the dominant terms 

for both limits, in the weak tension regime the fluctuation spectrum reads 

M = \f 4 - (25) 
crq^ + Kq^ 

Comparing this result with the fluctuation spectrum obtained for membranes with active proteins IT^ [TB^ , no 
novel nonequilibrium contribution to the height fluctuations is found here. The effect of the nonequilibrium reaction 
in stable membranes results in a change between a regime governed by Kg// to another one with the actual rigidity 
K. Thus, the reaction makes the membrane more rigid by simply removing the composition/curvature coupling effect 
that diminished the rigidity in the equilibrium situation. This change is more evident when the two components of 
the bilayer are rather different in shape and miscible but close to the critical temperature. In this situation (large 
and a < 0) the difference between n, and Kg/ / might be experimentally observable. 

Some caution, however, must be observed when deriving Eq. H25|l . In order to evaluate Eq. (|20|l for small 5, we 
considered the terms to be dominant with respect to the terms 2F {k.H'^ — a) q^ (subdominant) and also with 
respect to the subsubdominant quartic contributions (cr + [kHq — a)) q'^ (which are those that led to Eq. I|23|) '). 
However, the analysis becomes more difficult if one notices that we are not in the region of asymptotically small 
wavenumbers since the smallest q accessible to an experiment is of order 1/L, where L is the linear system size. Thus, 
one must look in detail at the values of the parameters in order to determine exactly which terms initially considered 
subdominant might indeed be dominant when the reaction is present. F^ is dominant if F > 2 [kHq — a) with 
qmin ~ ^/L and thus dependent on the system size. We know that for a sufficiently strong reactive process, F^ 
dominates and Eq. H25|l holds, but for intermediate values of F the spectrum could change significatively. When 
2F (kH§ — a) q^ becomes dominant, (|/igP) = 7 ^ — for a tensionless membrane. These regimes are 

captured in Fig.|51 where a set of curves is plotted for cr = 0. Very large and zero values for F show the q~'^ 

behavior, although with different rigidity factors {k and n^fj, respectively). At intermediate F, the q^'^ behavior 
appears as predicted above. 

The inclusion of hydrodynamic interactions can again be accomplished by considering A — {4:fiq) in the Fourier 
transform of Eq. H15() . This modifies the coefficients 021 and 022 in Eq. IjlSII . which are now divided by (4/15). 
Accordingly, the fh thermal noise correlations also change to {fh{j^, t)fh(j^' , t')) — ^^d{r-r')6{t-t'), and therefore 
the parameter A' in Eq. 1)22(1 has to be divided by (4/i(7) as well. However, with these modifications the evaluation of 
again leads to the results in Eqs. H23|) and H25|l. 



IV. CONCLUSIONS 



Starting with a simple model of a deformable reactive membrane composed of two differently shaped molecules, we 
show that stationary finite-sized patterns may appear under some parameter conditions for the immiscibility situation 
as a result of the competition between phase segregation and reaction. These structures involve heterogeneous distri- 
butions of composition and curvature whose sizes are determined by the nonequilibrium reactive process. For typical 
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FIG. 8: Log-log plot for evaluated from Eq. ^ in a tensionless membrane for different values of F. The other parameters 

are a — —0.01, k = 10 and Ho = 0.2 {neff = 2 in these cases). The slopes —1/4 and —1/3 are plotted to identify the different 
regimes. 



values of the viscosity of water and lipid lateral diffusion constants in bilayers, and at normal room temperatures, 
such patterns are predicted to have a size of a few microns (see the discussion at the end of Sec. Ill B|) . Therefore, 
this behavior would correspond to a reliable pattern formation mechanism in lipid membranes which we believe to 
be experimentally accessible in giant synthetic vesicles. The amplitude of these patterns is modulated by the bi- 
layer rigidity and the spontaneous curvature of its components. In our numerical calculations we have used realistic 
typical values for the rigidity, while the spontaneous curvature depends on the specific geometry of the membrane 
constituents. We specifically propose that azobenzene compounds, which are known to show amphiphilic behavior 
in Langmuir monolayers, and whose shapes are strongly modified by means of well-known photoisomerization reac- 
tions 38, 39], might be suitable to test our predictions. The selection of the applied light wavelength and its intensity 
may determine both the fraction of the two isomers (0o) and the strength of the reactive process (F), respectively. 
Thus, experimental work on synthetic vesicle membranes made of these compounds can be specifically designed to 
confirm the results of this model. 

In the same context, for miscible components membranes we have found a difference between the equilibrium 
and nonequilibrium situations. The effects of the composition/curvature coupling and the reactive process on the 
membrane rigidity are established. Micropipet experiments with the proposed membrane systems might confirm 
these results. 
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APPENDIX A: DERIVATION OF THE AMPLITUDE EQUATIONS 

The starting point of the analysis is the one dimensional version of the model for reactive membranes presented 
herein. 
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l^t = (Ki?o ^ Vxx + 3/3 + (Pxx + 6^ ((y5 + 0o) (<^2;)^ - l^xxxx 

(a -3/302)2 
47 

/it = —nhxxxx + nHQifxx- 

For convenience, we have introduced in Eqs. (|A1(I some notation changes and definitions with respect to the ones 
that appear in Eqs. Q). Thus, in Eqs. (|A1|I subscripts indicate partial derivatives, the field ip — (j) — (f>Q has been 
defined, and we have introduced the control parameter e = (Fc — F) /Fc that accounts for the "distance" to the 
bifurcation between a homogenous state (e < 0) and pattern formation (e > 0). The homogeneous state according to 
this definitions corresponds to ip — and arbitrary h — h. 

As shown in Sec. Ijll A|) . by linearizing Eqs. IjAip one can easily check that if e > then the homogeneous state 
becomes unstable and, 

(p{x) ^ A exp (iqcx) + A* exp {-iqcx) , ^^^^ 
h (x) — B exp (iqcx) + B* exp {—iqcx) , 

is a solution in the steady state if q"^ — (^a — 3/30g) / (27). It is worth noting that in Eqs. HA2|I we have arbitrarily 
taken h — without any loss of generality. Moreover, by substituting ljA2p into (|A1|I and expanding up to the first 
harmonic, O (exp (iqcx)), one finds that the amplitudes scale as a function of £ as A, _B ~ y^. Thus, we expect that 
near the bifurcation the following expansion holds, 

00 

(A3) 

n=l 

By computing the linear growth rate, LUq, that is, the largest eigenvalue of the linear problem (see Eqs. (0), we 
also note that, 

UJg\q^qa ~ Cie + C2 ((7 - qcf , 

where Ci are constants. Thus, as a function of e the width of the band of unstable modes scales as ~ e^/^^ Then, since 
all modes exp (iqx) can be written as exp (i (q — qc) x) exp {—iqcx), a separation of spatial scales can be performed 
between the most unstable mode (fast) and the rest of the modes of the unstable band (slow). Let us call the slow 
modulation spatial scale X, such that X = e^l'^x, where x will now stand for the fast spatial scale. We also note that 

exp {bJqt)\q-^qa ~ CXp {st) . 

Therefore we can define a slow time scale as a function of the control parameter, T = et. The separation of scales can 
be implemented in Eqs. IjAip by replacing the spatial and temporal derivatives according to the chain rule such that 
dx—' dx+ e^l'^dx and 9* edr- 

By implementing the separation of scales and substituting Eqs. (IA3() into Eqs. (|A1() we obtain a rather cumbersome 
expansion in terms of e. The lowest order contribution is of O (e^^^) and reads 

{ml + i?o« - «) - iv'^xLx - ("-y*^") yd) _ i/o^/iil. = 0, ^^^^ 

Note that Eqs. ljA4p correspond to the linearized version of Eqs. (|A1|I in the stationary state. We define the linear 
operator 

L = I ^^"^"^ + (3/3(/)§ + H^K - a) dxx - idxxxx -HqkOxxxx \ _ 

V HqkOxx ~l^dxxxx I 
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Then, Eqs. (|A4|) can be trivially written as hxi = 0, where (Xn)"^ = {^^"\ /i'"^). The contributions of the next order, 
0(e), are ]Lx2 = {{<f''^'^;h^^'^}) where = (4''\V'2^^' 



Finally, at order e^/^ we get Lxa = ips ({<p'^\ (/J^^^; ft.'^^}), where ^/il" = (^■03°\ V's 



reads, 



(a - H; 



3/3 



2.^ 



+ 27 (3^LW 



(2) 



2/l 



(2) 



+ 2/1^'^ 



+ 2^i^ 



We could continue up to any order with the expansion. In all cases we will obtain a nonlinear equation, such that 
at order e"/^, 



(A5) 



However, at order e^/^ we are already able to extract a closed evolution equation for the amplitudes of the pattern 
and so we will stop at that order. 

Our task is to solve the hierarchy of equations given by (jASp . At order e^/^ the problem is homogeneous and with 
appropriate boundary conditions, 



(^(1) = A {X, T) exp {iq^x) + A* {X, T) exp {-iqcx) , 
h^^^ = B (X, T) exp {iq^x) + B* {X, T) exp {-iq^x) , 



(A6) 



is a solution. However, the amplitudes A and B are undetermined at this point. The subsequent orders are no 
longer homogeneous and therefore their solvability can not be ensured unless one implements the so-called Fredholm 
alternative theorem In our case the application of the theorem simply states, as a recipe, that for Eqs. (|A5() 

to have a solution the functions tpn can not contain the fundamental mode exp{±iqcx). Thus, by substituting the 
solution (jA6p into the next order or the hierarchy and imposing the solvability condition we obtain 



3 (a - ml) 



(A {X, T) f exp ii2qcx) + {A (X, T)*f exp {-i2qcx) 
{A (X, T)f exp {i2q,x) + {A {X, T)*)' exp {-i2q,x) 



Once again, the value of A and B can not be determined at this order. However, at order e^/^ the application of 
the Fredholm theorem provides the conditions that determine the values of the amplitudes A, B. These conditions 
constitute the amplitude equations for our pattern forming system. 



drA 



47 

+ 2[a-3P^l + ^ 



27 

dxxA H ^ -dxxB, 



7 



(A7) 



3k (a - 3B6f,) 
drB = HoKdxxA + — ^ ^-^dxxB. 
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Equations (|A7|) can be rewritten in terms of x and t to readily obtain Eqs. 
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